Import and get ready:

Import previous metaphlan data so we will get the metadata from there:

load(here::here("../../data/processed_data/metaphlan/01_data.Rdata"))

Define path of the PWY and Gene family table.- output of HumAn

pathway <- "/Users/ljc444/Documents/Caroline/temp_results/humann/Merged_pathabundance_cpm.tsv"
kegg <- "/Users/ljc444/Documents/Caroline/temp_results/humann/Merged_genefamilies_uniref90_ko_renamed_relab.tsv"
l4c <- "/Users/ljc444/Documents/Caroline/temp_results/humann/Merged_genefamilies_uniref90_level4ec_renamed_relab.tsv"
go <- "/Users/ljc444/Documents/Caroline/temp_results/humann/Merged_genefamilies_uniref90_go_renamed_relab.tsv"
eggnog <- "/Users/ljc444/Documents/Caroline/temp_results/humann/Merged_genefamilies_uniref90_eggnog_relab.tsv"
rxn <- "/Users/ljc444/Documents/Caroline/temp_results/humann/Merged_genefamilies_uniref90_rxn_renamed_relab.tsv"
pfam <- "/Users/ljc444/Documents/Caroline/temp_results/humann/Merged_genefamilies_uniref90_pfam_renamed_relab.tsv"

Arguments of the import_humann_pathway_mia function:

#' Import and Process HUMAnN Pathway Data for Microbial Ecology Analysis
#'
#' This function imports HUMAnN pathway data, processes it into a `phyloseq` object, 
#' and cleans and annotates the data to facilitate downstream analysis in microbial ecology studies.
#' 
#' @param pathway A file path or object containing HUMAnN pathway or gene family level data.
#' @param rank The taxonomic rank used for filtering or organizing the taxonomy table. Default is `"pathway"`.
#' @param tax_clean_gsub A regular expression pattern to clean unwanted characters or prefixes in taxonomy labels. Default is `"[a-t]__"`.
#' @param no_sel A vector of values to exclude during filtering, such as `NA` or `"unclassified"`. Default is `c(NA, "<NA>", "NA", "unclassified")`.
#' @param sub_pat_a A string pattern for initial cleaning of sample names. Default is `"_trimmed_HF_merged_Abundance-RELAB"`.
#' @param sub_pat_b A string pattern for additional cleaning of sample names. Default is `"_subset50"`.
#' @param tax_na_if A string indicating taxonomy entries to replace with `NA`. Default is `"unclassified"`.
#' @param sample_name_clean_a A regular expression pattern for extracting meaningful parts of sample names. Default is `"[^_]+"`.
#' @param sample_name_clean_b A regular expression for removing specific patterns in sample names. Default is `"[A-Z]"`.
#' @param meta File path to metadata (in Excel format). Default points to a metadata file for the Deerland Study.
#' @param sample_column The column name in the metadata file corresponding to sample IDs. Default is `"ID"`.
#' @param sample_nam_prefix A prefix to add to sample names. Default is `"S_"`.
#' @param rm_un Logical; whether to remove taxonomy entries starting with `"UN"`. Default is `TRUE`.
#' @param return_unstrat Logical; whether to return only unstratified pathways. Default is `TRUE`.
#' @param str_trun An integer specifying the maximum length of pathway names for truncation. Default is `40`.
#' @param add_MetaCyc_pathway_map Logical; whether to include MetaCyc pathway mapping. Default is `TRUE`.
#' @param add_CHOCOPhlAn_taxonomy Logical; whether to annotate with CHOCOPhlAn taxonomy. Default is `TRUE`.
#' @param final_tax_table_ranks A vector of column names defining the final structure of the taxonomy table. Default includes taxonomic and pathway features.
#' 
#' @return A `phyloseq` object with processed HUMAnN pathway data, cleaned sample names, taxonomy table, and metadata.
#' 
#' @details 
#' This function performs the following operations:
#' - Imports HUMAnN pathway data using `mia::importHUMAnN`.
#' - Converts the data into a `phyloseq` object for downstream microbial ecology analysis.
#' - Cleans sample names based on specified substitution patterns and regular expressions.
#' - Adds metadata to the `phyloseq` object for sample-level annotations.
#' - Processes and cleans the taxonomy table, replacing "unclassified" entries with `NA`, and optionally removing entries starting with "UN".
#' - Annotates pathways using MetaCyc and/or CHOCOPhlAn taxonomy data, with abbreviations for long pathway names and features.
#' 
#' @examples
#' \dontrun{
#' pathway_data <- "path/to/humann_pathway_data.tsv"
#' processed_data <- import_humann_pathway_mia(
#'   pathway = pathway_data,
#'   meta = "path/to/metadata.xlsx",
#'   sample_column = "SampleID"
#' )
#' }
#' 

Extract the metadata from our metaphlan phyloseq object:

ps_up %>% 
  sample_data() %>% 
  data.frame() %>% 
  rownames_to_column('tmp') -> metaph_metadata

Import HumAnN partway data and combine with the metadata of our metaphlan object:

pathway %>% 
  import_humann_pathway_mia(meta = metaph_metadata, 
                            sample_column = 'tmp',
                            rm_un = FALSE, # Remove UNMAPPED and UNINTEGRATED Pathways
                            return_unstrat = FALSE, # 
                            add_CHOCOPhlAn_taxonomy = TRUE, # add Metaphlan/humann taxonomic path 
                            add_MetaCyc_pathway_map = FALSE) -> pw_mia # add MetaCyc pathway hierarchy

pw_mia %>% 
  tax_table() %>% 
  head()
## Taxonomy Table:     [ 6 taxa by 8 taxonomic ranks ]:
##                          feature Kingdom Phylum Class Order Family Genus Species
##                          <chr>   <chr>   <chr>  <chr> <chr> <chr>  <chr> <chr>  
## 1 UNMAPPED               UNMAPP… <NA>    <NA>   <NA>  <NA>  <NA>   <NA>  <NA>   
## 2 UNINTEGRATED           UNINTE… <NA>    <NA>   <NA>  <NA>  <NA>   <NA>  <NA>   
## 3 UNINTEGRATED|g__Abiot… UNINTE… Bacter… Firmi… Baci… Lact… Aeroc… Abio… Abiotr…
## 4 UNINTEGRATED|g__Abiot… UNINTE… Bacter… Firmi… Baci… Lact… Aeroc… Abio… Abiotr…
## 5 UNINTEGRATED|g__Actin… UNINTE… Bacter… Actin… Acti… Acti… Actin… Acti… Actino…
## 6 UNINTEGRATED|g__Actin… UNINTE… Bacter… Actin… Acti… Acti… Actin… Acti… Actino…
pw_mia %>% 
  tax_table() %>% 
  tail()
## Taxonomy Table:     [ 6 taxa by 8 taxonomic ranks ]:
##                          feature Kingdom Phylum Class Order Family Genus Species
##                          <chr>   <chr>   <chr>  <chr> <chr> <chr>  <chr> <chr>  
## 1 VALSYN-PWY: L-valine … VALSYN… Bacter… Firmi… Nega… Veil… Veill… Veil… Veillo…
## 2 VALSYN-PWY: L-valine … VALSYN… Bacter… Firmi… Nega… Veil… Veill… Veil… Veillo…
## 3 VALSYN-PWY: L-valine … VALSYN… Bacter… Firmi… Nega… Veil… Veill… Veil… Veillo…
## 4 VALSYN-PWY: L-valine … VALSYN… Bacter… Firmi… Nega… Veil… Veill… Veil… Veillo…
## 5 VALSYN-PWY: L-valine … VALSYN… Bacter… Firmi… Nega… Veil… Veill… Veil… Veillo…
## 6 VALSYN-PWY: L-valine … VALSYN… <NA>    <NA>   <NA>  <NA>  <NA>   <NA>  <NA>

By default return Unstratified data:

pathway %>% 
  import_humann_pathway_mia(meta = metaph_metadata, 
                            sample_column = 'tmp') -> pw_mia

pw_mia %>% 
  tax_table() %>% 
  head()
## Taxonomy Table:     [ 6 taxa by 10 taxonomic ranks ]:
##                Superclass1 Superclass2 feature Kingdom Phylum Class Order Family
##                <chr>       <chr>       <chr>   <chr>   <chr>  <chr> <chr> <chr> 
## 1 1CMET2-PWY   Bios.       Cofactor, … folate… <NA>    <NA>   <NA>  <NA>  <NA>  
## 2 ALLANTOINDE… Degrad./Ut… Amide, Ami… SPWY. … <NA>    <NA>   <NA>  <NA>  <NA>  
## 3 ANAEROFRUCA… Gen. of Pr… Ferm.       homola… <NA>    <NA>   <NA>  <NA>  <NA>  
## 4 ANAGLYCOLYS… Gen. of Pr… Glycol.     Glycol… <NA>    <NA>   <NA>  <NA>  <NA>  
## 5 ARG+POLYAMI… Bios.       Amide, Ami… SPWY. … <NA>    <NA>   <NA>  <NA>  <NA>  
## 6 ARGININE-SY… Bios.       Amino Acid… L-orni… <NA>    <NA>   <NA>  <NA>  <NA>  
## # ℹ 2 more taxonomic ranks: Genus <chr>, Species <chr>

We can also import Gene family data e.g., KEGG:

kegg %>% 
  import_humann_pathway_mia(meta = metaph_metadata, 
                            sample_column = 'tmp',
                            rm_un = TRUE, 
                            return_unstrat = FALSE,
                            add_CHOCOPhlAn_taxonomy = TRUE, 
                            add_MetaCyc_pathway_map = FALSE) -> kegg_mia # add_MetaCyc_pathway_map only for pathway level data

kegg_mia %>% 
  tax_table() %>% 
  head()
## Taxonomy Table:     [ 6 taxa by 8 taxonomic ranks ]:
##                          feature Kingdom Phylum Class Order Family Genus Species
##                          <chr>   <chr>   <chr>  <chr> <chr> <chr>  <chr> <chr>  
## 1 K00001: alcohol dehyd… K00001… <NA>    <NA>   <NA>  <NA>  <NA>   <NA>  <NA>   
## 2 K00001: alcohol dehyd… K00001… Bacter… Prote… Beta… Burk… Comam… Delf… Delfti…
## 3 K00001: alcohol dehyd… K00001… Bacter… Prote… Beta… Burk… Comam… Delf… Delfti…
## 4 K00001: alcohol dehyd… K00001… Bacter… Prote… Beta… Burk… Comam… Delf… Delfti…
## 5 K00001: alcohol dehyd… K00001… Bacter… Firmi… Baci… Lact… Lacto… Lact… Lactob…
## 6 K00001: alcohol dehyd… K00001… Bacter… Firmi… Baci… Lact… Lacto… Lact… Lactob…

Let’s focus on the unstratified data, then we can always go back to the the stratified data to check which taxa are contributing to particular functions:

kegg %>% 
  import_humann_pathway_mia(meta = metaph_metadata, 
                            sample_column = 'tmp',
                            rm_un = TRUE, 
                            return_unstrat = TRUE,
                            add_CHOCOPhlAn_taxonomy = TRUE, 
                            add_MetaCyc_pathway_map = FALSE) -> kegg_mia

kegg_mia %>% 
  tax_table() %>% 
  head()
## Taxonomy Table:     [ 6 taxa by 8 taxonomic ranks ]:
##          feature                 Kingdom Phylum Class Order Family Genus Species
##          <chr>                   <chr>   <chr>  <chr> <chr> <chr>  <chr> <chr>  
## 1 K00001 K00001: alcohol dehydr… <NA>    <NA>   <NA>  <NA>  <NA>   <NA>  <NA>   
## 2 K00002 K00002: alcohol dehydr… <NA>    <NA>   <NA>  <NA>  <NA>   <NA>  <NA>   
## 3 K00003 K00003: homoserine deh… <NA>    <NA>   <NA>  <NA>  <NA>   <NA>  <NA>   
## 4 K00004 K00004: (R,R)-butanedi… <NA>    <NA>   <NA>  <NA>  <NA>   <NA>  <NA>   
## 5 K00005 K00005: glycerol dehyd… <NA>    <NA>   <NA>  <NA>  <NA>   <NA>  <NA>   
## 6 K00006 K00006: glycerol-3-pho… <NA>    <NA>   <NA>  <NA>  <NA>   <NA>  <NA>

Use previous function for alpha - beta - differentially abundant taxa/features (e.g., Pathway, gene) tests:

Alpha:

PWY:

Let’s focus on the unstratified data, then we can always go back to the the stratified data to check which taxa are contributing to particular functions:

alpha_measures = c("observed", "diversity_shannon", "diversity_inverse_simpson", "evenness_pielou")

pw_mia %>% 
  transform_sample_counts(function(x) x/sum(x) * 100) %>% 
  phyloseq_alphas() %>% 
  pivot_longer(cols = all_of(alpha_measures), values_to = "value", names_to = 'alphadiversiy', values_drop_na  = TRUE) %>%
  mutate(alphadiversiy = fct_relevel(alphadiversiy, alpha_measures)) %>% 
  phyloseq_explore_alpha(facet_formula = "alphadiversiy ~ Sample",
                         group_by_stats = c("alphadiversiy", "Sample"),
                         stat_formula = "value ~ Time",
                         padjust_method = "fdr") -> out_pway

ls(out_pway)
## [1] "alpha_legend"  "alpha_plot"    "alpha_summary" "anova_result" 
## [5] "anova_table"   "kruskal_test"  "t_test"        "wilcox_test"

Metrics are not calculated on the taxa here but on the functional Pathways:

out_pway$alpha_plot

KO:

We can do the same with Gene level information, e.g.: Kegg’s Orthologues:

kegg_mia %>%
  physeq_glom_rename(taxrank = "feature", 
                     rename_ASV = "feature") %>% 
  transform_sample_counts(function(x) x/sum(x) * 100) %>% 
  phyloseq_alphas() %>% 
  pivot_longer(cols = all_of(alpha_measures), values_to = "value", names_to = 'alphadiversiy', values_drop_na  = TRUE) %>%
  mutate(alphadiversiy = fct_relevel(alphadiversiy, alpha_measures)) %>% 
  phyloseq_explore_alpha(facet_formula = "alphadiversiy ~ Sample",
                         group_by_stats = c("alphadiversiy", "Sample"),
                         stat_formula = "value ~ Time",
                         padjust_method = "fdr",
                         # stat_paired = FALSE,
                         ref_group_stat = "TP1") -> out_kegg

out_kegg$alpha_plot

Beta:

PWY:

out$PCOA

out$pcoas

out$Sample$Saliva

out$dist_box

KO:

kegg_mia %>%
  physeq_glom_rename(taxrank = "feature", 
                     rename_ASV = "feature")  -> ps_tmp

ps_tmp %>%
  # subset_taxa(Class != "UNCLASSIFIED") %>% 
  transform_sample_counts(function(x) x/sum(x) * 100) %>%
  phyloseq_compute_bdiv() -> beta


ps_tmp %>% 
  # subset_taxa(Class != "UNCLASSIFIED") %>% 
  transform_sample_counts(function(x) x/sum(x) * 100) %>%
  microViz::dist_calc(., dist = "robust.aitchison") %>% 
  microViz::dist_get() %>% 
  magrittr::divide_by(100) -> beta$rAitchison

beta$bray <- NULL
beta$sorensen <- NULL

ps_tmp %>%
  # subset_taxa(Class != "UNCLASSIFIED") %>% 
  transform_sample_counts(function(x) x/sum(x) * 100) %>%
  phyloseq_explore_beta(ps_up = .,
                        beta = beta,
                        color_group = "Time",
                        shape_group = "Sample",
                        alpha = NULL,
                        col_pal = time_pal,
                        fill_pal = time_pal,
                        path_group = "interaction(Sample,Subject)",
                        facet_formula = "Sample ~ .",
                        axis1 = 1,
                        axis2 = 2,
                        seed = 123,
                        permanova_terms = c("Sample", "Time"),
                        strata = "none", 
                        metadata_dist_boxplot = c("Subject", "Time", "Sample")) -> out
out$Sample
## $Saliva

## 
## $Plaque

Feature visualisation/ differential abundance tests:

Since some of the functions for differential abundance analyses and taxa/feature visualization expect to find Kingdom, Phylum, Class, Order, Family, Genus and Species as rank_names we can check the structure of the data and the name of the columns organizing that hierarchy:

pw_mia %>% 
  tax_table() %>% 
  head()
## Taxonomy Table:     [ 6 taxa by 10 taxonomic ranks ]:
##                Superclass1 Superclass2 feature Kingdom Phylum Class Order Family
##                <chr>       <chr>       <chr>   <chr>   <chr>  <chr> <chr> <chr> 
## 1 1CMET2-PWY   Bios.       Cofactor, … folate… <NA>    <NA>   <NA>  <NA>  <NA>  
## 2 ALLANTOINDE… Degrad./Ut… Amide, Ami… SPWY. … <NA>    <NA>   <NA>  <NA>  <NA>  
## 3 ANAEROFRUCA… Gen. of Pr… Ferm.       homola… <NA>    <NA>   <NA>  <NA>  <NA>  
## 4 ANAGLYCOLYS… Gen. of Pr… Glycol.     Glycol… <NA>    <NA>   <NA>  <NA>  <NA>  
## 5 ARG+POLYAMI… Bios.       Amide, Ami… SPWY. … <NA>    <NA>   <NA>  <NA>  <NA>  
## 6 ARGININE-SY… Bios.       Amino Acid… L-orni… <NA>    <NA>   <NA>  <NA>  <NA>  
## # ℹ 2 more taxonomic ranks: Genus <chr>, Species <chr>

Here ‘Superclass1’, ‘Superclass2’ and ‘feature’ are the interesting ones:

So we can use tax_mutate to organise (and copy) those data into Kingdom, Phylum, Class, Order, Family, Genus and Species columns.

pw_mia %>% 
  tax_mutate(Kingdom = "MetaCyc", 
             Phylum = Superclass1,
             Class = Superclass1,
             Order = Superclass1,
             Family = Superclass2,
             Genus = Superclass2,
             Species = feature) %>% 
  physeq_sel_tax_table(c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")) -> ps_tmp

Then the next step looks familiar!

ps_tmp %>% 
  phyloseq_top_heatmap_barplot(facet_formula = "Sample_Type ~ Time" , group_var = "Sample_Time",
                               tax_levels = c("Phylum", "Genus", "Species"),
                               ntax = 4, ntax_species = 20, plot_heights = c(2.2, 1.5, 6),
                               boxplot_main_group = "Genus",barplot_level = "Species",
                               rm_unclassified = FALSE,
                               facet_by = c("Sample_Type", "Time"),
                               group_by = c("Sample_Type", "Time"),
                               facet_heat = "~ Sample_Type + Time ") -> t

t$heat_all

t$p

t$nested_legend

PWY:

ps_tmp %>% 
  subset_samples(Time == "TP1") %>% 
  phyloseq_diff_abundance(ps_tmp = ., 
                          approach = c("run_lefse"),
                          comp_group = "Sample_Type",
                          glom = "Species", # agglomerate data at the species level -> pathway
                          taxa_rank = "Species", # perform comparison at the Species level -> pathway
                          lefse_lda_cutoff = 3, # filter results so only strong differences are displayed 
                          pvalue_cutoff = 0.0001, # filter results so only strong differences are displayed 
                          palette = sample_pal) -> lefse_TP1
## Loading required package: microbiomeMarker
## Registered S3 methods overwritten by 'proxy':
##   method               from    
##   print.registry_field registry
##   print.registry_entry registry
## Registered S3 method overwritten by 'gplots':
##   method         from     
##   reorder.factor DescTools
## 
## Attaching package: 'microbiomeMarker'
## The following objects are masked from 'package:microbiome':
## 
##     abundances, aggregate_taxa
## The following object is masked from 'package:speedyseq':
## 
##     plot_heatmap
## The following object is masked from 'package:phyloseq':
## 
##     plot_heatmap
## Loading required package: magrittr
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:purrr':
## 
##     set_names
## The following object is masked from 'package:tidyr':
## 
##     extract
## Loading required package: ANCOMBC
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
## Welcome to compositions, a package for compositional data analysis.
## Find an intro with "? compositions"
## 
## Attaching package: 'compositions'
## The following object is masked from 'package:microbiomeMarker':
## 
##     normalize
## The following object is masked from 'package:ape':
## 
##     balance
## The following objects are masked from 'package:stats':
## 
##     anova, cor, cov, dist, var
## The following object is masked from 'package:graphics':
## 
##     segments
## The following objects are masked from 'package:base':
## 
##     %*%, norm, scale, scale.default
lefse_TP1$mmlefse_p

Gene Families:

kegg_mia %>% 
  tax_mutate(Kingdom = "KEGG", 
             Phylum = feature,
             Class = feature,
             Order = feature,
             Family = feature,
             Genus = feature,
             Species = feature) %>% 
  physeq_sel_tax_table(c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")) -> ps_tmp

Then the next step looks familiar!

ps_tmp %>% 
  phyloseq_top_heatmap_barplot(facet_formula = "Sample_Type ~ Time" , group_var = "Sample_Time",
                               tax_levels = c("Phylum", "Genus", "Species"),
                               ntax = 4, ntax_species = 22, plot_heights = c(2.2, 1.5, 6),
                               boxplot_main_group = "Family",
                               rm_unclassified = FALSE,
                               facet_by = c("Sample_Type", "Time"),
                               group_by = c("Sample_Type", "Time"),
                               facet_heat = "~ Sample_Type + Time ") -> t

t$heat$Species

ps_tmp %>% 
  subset_samples(Time == "TP1") %>% 
  phyloseq_diff_abundance(ps_tmp = ., 
                          approach = c("run_lefse"),
                          comp_group = "Sample_Type",
                          glom = "Species", # agglomerate data at the species level -> pathway
                          taxa_rank = "Species", # perform comparison at the Species level -> pathway
                          lefse_lda_cutoff = 3, # filter results so only strong differences are displayed 
                          pvalue_cutoff = 0.0001, # filter results so only strong differences are displayed 
                          palette = sample_pal) -> lefse_TP1


lefse_TP1$mmlefse_p

ps_tmp %>% 
  subset_samples(Sample_Type == "Plaque" & Time%in% c("TP1", "TP2")) %>% 
  phyloseq_diff_abundance(ps_tmp = ., 
                          approach = c("run_lefse"),
                          comp_group = "Time",
                          glom = "Species", # agglomerate data at the species level -> pathway
                          taxa_rank = "Species", # perform comparison at the Species level -> pathway
                          lefse_lda_cutoff = 2.5, # filter results so only strong differences are displayed 
                          pvalue_cutoff = 0.001, # filter results so only strong differences are displayed 
                          palette = time_pal) -> lefse_TP1


lefse_TP1$mmlefse_p

sessionInfo()
## R version 4.4.0 (2024-04-24)
## Platform: aarch64-apple-darwin20
## Running under: macOS 15.1.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Europe/Zurich
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] compositions_2.0-8      nlme_3.1-166            ANCOMBC_2.8.0          
##  [4] magrittr_2.0.3          microbiomeMarker_1.10.0 ggnested_0.1.1         
##  [7] vegan_2.6-8             lattice_0.22-6          permute_0.9-7          
## [10] GUniFrac_1.8            ape_5.8                 microbiome_1.28.0      
## [13] metagMisc_0.5.0         rstatix_0.7.2           microViz_0.12.4        
## [16] ampvis2_2.8.9           ggpubr_0.6.0            file2meco_0.9.0        
## [19] microeco_1.10.1         reshape2_1.4.4          scales_1.3.0           
## [22] speedyseq_0.5.3.9021    phyloseq_1.50.0         readxl_1.4.3           
## [25] lubridate_1.9.3         forcats_1.0.0           stringr_1.5.1          
## [28] dplyr_1.1.4             purrr_1.0.2             readr_2.1.5            
## [31] tidyr_1.3.1             tibble_3.2.1            ggplot2_3.5.1          
## [34] tidyverse_2.0.0        
## 
## loaded via a namespace (and not attached):
##   [1] coin_1.4-3                      IRanges_2.39.2                 
##   [3] gld_2.6.6                       nnet_7.3-19                    
##   [5] Biostrings_2.73.2               TH.data_1.1-2                  
##   [7] vctrs_0.6.5                     energy_1.7-12                  
##   [9] digest_0.6.37                   png_0.1-8                      
##  [11] shape_1.4.6.1                   proxy_0.4-27                   
##  [13] Exact_3.3                       registry_0.5-1                 
##  [15] rbiom_1.0.3                     ggrepel_0.9.6                  
##  [17] mediation_4.5.0                 plotwidgets_0.5.1              
##  [19] MASS_7.3-61                     foreach_1.5.2                  
##  [21] BiocGenerics_0.52.0             withr_3.0.2                    
##  [23] xfun_0.49                       ggfun_0.1.7                    
##  [25] survival_3.7-0                  doRNG_1.8.6                    
##  [27] ggbeeswarm_0.7.2                gmp_0.7-5                      
##  [29] tidytree_0.4.6                  zoo_1.8-12                     
##  [31] GlobalOptions_0.1.2             gtools_3.9.5                   
##  [33] DEoptimR_1.1-3                  Formula_1.2-5                  
##  [35] httr_1.4.7                      rhdf5filters_1.17.0            
##  [37] rhdf5_2.49.0                    rstudioapi_0.17.1              
##  [39] UCSC.utils_1.2.0                generics_0.1.3                 
##  [41] base64enc_0.1-3                 S4Vectors_0.43.2               
##  [43] zlibbioc_1.51.2                 ScaledMatrix_1.14.0            
##  [45] ca_0.71.1                       statip_0.2.3                   
##  [47] GenomeInfoDbData_1.2.13         SparseArray_1.5.45             
##  [49] ade4_1.7-22                     doParallel_1.0.17              
##  [51] evaluate_1.0.1                  S4Arrays_1.5.11                
##  [53] glmnet_4.1-8                    hms_1.1.3                      
##  [55] GenomicRanges_1.57.2            irlba_2.3.5.1                  
##  [57] colorspace_2.1-1                modeltools_0.2-23              
##  [59] viridis_0.6.5                   ggtree_3.14.0                  
##  [61] robustbase_0.99-4-1             DECIPHER_3.1.5                 
##  [63] scuttle_1.15.5                  cowplot_1.1.3                  
##  [65] matrixStats_1.4.1               class_7.3-22                   
##  [67] Hmisc_5.2-0                     pillar_1.9.0                   
##  [69] iterators_1.0.14                decontam_1.26.0                
##  [71] plotROC_2.3.1                   caTools_1.18.3                 
##  [73] compiler_4.4.0                  beachmat_2.22.0                
##  [75] stringi_1.8.4                   biomformat_1.34.0              
##  [77] DescTools_0.99.57               TSP_1.2-4                      
##  [79] stabledist_0.7-2                minqa_1.2.8                    
##  [81] SummarizedExperiment_1.36.0     plyr_1.8.9                     
##  [83] crayon_1.5.3                    abind_1.4-8                    
##  [85] scater_1.34.0                   timeSeries_4041.111            
##  [87] gridGraphics_0.5-1              ggtext_0.1.2                   
##  [89] locfit_1.5-9.10                 bit_4.5.0                      
##  [91] mia_1.14.0                      rootSolve_1.8.2.4              
##  [93] sandwich_3.1-1                  libcoin_1.0-10                 
##  [95] codetools_0.2-20                multcomp_1.4-26                
##  [97] BiocSingular_1.21.4             bslib_0.8.0                    
##  [99] slam_0.1-54                     e1071_1.7-16                   
## [101] lmom_3.2                        GetoptLong_1.0.5               
## [103] plotly_4.10.4                   multtest_2.61.0                
## [105] MultiAssayExperiment_1.32.0     metagenomeSeq_1.46.0           
## [107] splines_4.4.0                   circlize_0.4.16                
## [109] Rcpp_1.0.13-1                   sparseMatrixStats_1.17.2       
## [111] cellranger_1.1.0                gridtext_0.1.5                 
## [113] knitr_1.48                      utf8_1.2.4                     
## [115] here_1.0.1                      clue_0.3-65                    
## [117] lme4_1.1-35.5                   fBasics_4041.97                
## [119] fs_1.6.5                        checkmate_2.3.2                
## [121] DelayedMatrixStats_1.28.0       Rdpack_2.6.1                   
## [123] expm_1.0-0                      ggplotify_0.1.2                
## [125] gsl_2.1-8                       ggsignif_0.6.4                 
## [127] Matrix_1.7-1                    statmod_1.5.0                  
## [129] tzdb_0.4.0                      lpSolve_5.6.21                 
## [131] pkgconfig_2.0.3                 tools_4.4.0                    
## [133] cachem_1.1.0                    rbibutils_2.3                  
## [135] viridisLite_0.4.2               DBI_1.2.3                      
## [137] numDeriv_2016.8-1.1             rmutil_1.1.10                  
## [139] fastmap_1.2.0                   rmarkdown_2.29                 
## [141] grid_4.4.0                      broom_1.0.7                    
## [143] sass_0.4.9                      patchwork_1.3.0                
## [145] stable_1.1.6                    carData_3.0-5                  
## [147] rpart_4.1.23                    farver_2.1.2                   
## [149] mgcv_1.9-1                      yaml_2.3.10                    
## [151] bayesm_3.1-6                    MatrixGenerics_1.18.0          
## [153] spatial_7.3-17                  foreign_0.8-87                 
## [155] cli_3.6.3                       stats4_4.4.0                   
## [157] lifecycle_1.0.4                 Biobase_2.66.0                 
## [159] mvtnorm_1.3-2                   bluster_1.15.1                 
## [161] backports_1.5.0                 modeest_2.4.0                  
## [163] BiocParallel_1.39.0             timechange_0.3.0               
## [165] gtable_0.3.6                    rjson_0.2.23                   
## [167] limma_3.61.12                   parallel_4.4.0                 
## [169] CVXR_1.0-14                     jsonlite_1.8.9                 
## [171] bitops_1.0-9                    seriation_1.5.6                
## [173] bit64_4.5.2                     Rtsne_0.17                     
## [175] yulab.utils_0.1.7               BiocNeighbors_1.99.3           
## [177] TreeSummarizedExperiment_2.14.0 RcppParallel_5.1.9             
## [179] jquerylib_0.1.4                 highr_0.11                     
## [181] timeDate_4041.110               lazyeval_0.2.2                 
## [183] htmltools_0.5.8.1               glue_1.8.0                     
## [185] fantaxtic_0.2.1                 Wrench_1.24.0                  
## [187] XVector_0.45.0                  rprojroot_2.0.4                
## [189] treeio_1.30.0                   gridExtra_2.3                  
## [191] boot_1.3-31                     igraph_2.1.1                   
## [193] R6_2.5.1                        gplots_3.2.0                   
## [195] DESeq2_1.45.3                   SingleCellExperiment_1.28.0    
## [197] Rmpfr_0.9-5                     labeling_0.4.3                 
## [199] cluster_2.1.6                   rngtools_1.5.2                 
## [201] Rhdf5lib_1.27.0                 aplot_0.2.3                    
## [203] GenomeInfoDb_1.42.0             nloptr_2.1.1                   
## [205] DirichletMultinomial_1.47.2     DelayedArray_0.31.14           
## [207] tidyselect_1.2.1                vipor_0.4.7                    
## [209] tensorA_0.36.2.1                htmlTable_2.4.3                
## [211] xml2_1.3.6                      inline_0.3.19                  
## [213] car_3.1-3                       KernSmooth_2.23-24             
## [215] rsvd_1.0.5                      munsell_0.5.1                  
## [217] data.table_1.16.2               htmlwidgets_1.6.4              
## [219] ComplexHeatmap_2.22.0           RColorBrewer_1.1-3             
## [221] rlang_1.1.4                     lmerTest_3.1-3                 
## [223] fansi_1.0.6                     beeswarm_0.4.0